#ifndef _MATRIX_H
#define _MATRIX_H
#include<bits/stdc++.h>
using namespace std;

template <class T> class matrix;
template <class T> class rowvec;
template <class T> class colvec;

template <class T> T operator * (const rowvec<T>&,const colvec<T>&);

template <class T> class matrix{
private:
	int n, m;
	T *a;
public:
	matrix(){
		a = NULL, n = m = 0;
	}
	matrix(int _n,int _m){
		if(_n <= 0 || _m <= 0) throw "Illegal Creating!";
		n = _n, m = _m;
		a = new T[_n * _m];
		memset(a, 0, sizeof(T) * _n * _m);
	}
	inline const T& at(const int& i)const{return a[i];}
	inline T& at(const int& i){return a[i];}
	inline const T* operator [] (const int& i)const{
		if(i < 0 || i >= n) throw "Illegal Visiting!";
		return &a[i * m];
	}
	inline T* operator [] (const int& i){
		if(i < 0 || i >= n) throw "Illegal Visiting!";
		return &a[i * m];
	}
	inline pair <int,int> size()const{
		return make_pair(n, m);
	}
	inline bool issq()const{
		return n == m;
	}
	inline void read(){
		for(int i = 0; i < n * m; ++i) cin >> a[i];
	}
	inline void write()const{
		for(int i = 0; i < n; ++i, cout << endl)
			for(int j = 0; j < m; ++j)
				cout << a[i * m + j] << " ";
	}
	matrix(const matrix& p){
		n = p.n, m = p.m;
		if(n && m) a = new T[n * m];
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < m; ++j)
				a[i * m + j] = p[i][j];
	}
	const matrix& operator = (const matrix& p){
		n = p.n, m = p.m;
		if(a != NULL) delete a;
		if(n && m) a = new T[n * m];
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < m; ++j)
				a[i * m + j] = p[i][j];
		return *this;
	}
	inline ~matrix(){
		if(n && m && a != NULL) delete a;
	}
	friend matrix<T> operator + (const matrix<T>& a,const matrix<T>& b){
		if(a.n != b.n || a.m != b.m)
			throw "Error: Illegal Addition!";
		int n = a.n, m = a.m;
		matrix<T> c = matrix <T> (n,m);
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < m; ++j)
				c[i][j] = a[i][j] + b[i][j];
		return c;
	}
	friend matrix<T> operator - (const matrix<T>& a, const matrix<T>& b){
		if(a.n != b.n || a.m != b.m)
			throw "Error: Illegal Subtract!";
		int n = a.n, m = a.m;
		matrix<T> c = matrix<T> (n,m);
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < m;++j)
				c[i][j] = a[i][j] - b[i][j];
		return c;
	}
	friend matrix<T> operator * (const matrix<T>& a, const matrix<T>& b){
		if(a.m != b.n)
			throw "Error: Illegal Multiply!";
		int n = a.n, p = a.m, m = b.m;
		matrix<T> c = matrix <T> (n,m);
		for(int k = 0; k < p; ++k)
			for(int i = 0; i < n; ++i)
				for(int j = 0; j < m; ++j)
					c[i][j] += a[i][k] * b[k][j];
		return c;
	}
	#define eps 1e-12
	inline T det()const{
		if(!issq())
			throw "Error: Non-Square Matrix have no determinent!";
		matrix<T> b = *this;
		T res = 1;
		for(int i = 0; i < n; ++i){
			int j = i;
			while(j < n && fabs(b[j][i]) < eps) ++j;
			if(j == n)return 0;
			if(j > i){
				res = -res;
				for(int k = i; k < n; ++k)
					swap(b[i][k], b[j][k]);
			}
			res *= b[i][i];
			for(int j = i + 1; j < n; ++j){
				T tmp = b[j][i] / b[i][i];
				for(int k = i; k < n; ++k)
					b[j][k] -= b[i][k] * tmp;
			}
		}
		return res;
	}
	inline matrix<T> inv()const{
		if(!issq())
			throw "Error: Non-Square Matrix have no Inverse!";
		if(fabs(det())<eps)
			throw "Error: Singular Matrix have no Inverse!";
		matrix<T> b = matrix<T> (n, n * 2);
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < n; ++j)
				b[i][j] = a[i*m+j];
		for(int i = 0; i < n; ++i)
			b[i][i+n] = 1;
		for(int i = 0; i < n; ++i){
			int j = i;
			while(j < n && fabs(b[j][i]) < eps) ++j;
			if(j > i){
				for(int k = i; k < n * 2; ++k)
					swap(b[i][k], b[j][k]);
			}
			T tmp = b[i][i];
			for(int j = 0; j < n; ++j)if(i ^ j){
				T z = b[j][i] / b[i][i];
				for(int k = 0; k < n * 2; ++k)
					b[j][k] -= b[i][k] * z;
			}
			for(int j = 0; j < n * 2;++j)
				b[i][j] /= tmp;
		}
		matrix<T> res = matrix<T> (n, n);
		for(int i = 0; i < n; ++i)
			for(int j = 0; j < n; ++j)
				res[i][j] = b[i][j+n];
		return res;
	}
};

template <class T> class rowvec : public matrix<T>{
public:
	inline rowvec():matrix<T>(){}
	inline rowvec(int _n):matrix<T>(1, _n){}
	inline const T operator [] (const int& i)const{
		return matrix<T>::at(i);
	}
	inline T& operator [] (const int& i){
		return matrix<T>::at(i);
	}
	inline int size()const{
		return matrix<T>::size().second;
	}
	friend T operator * (const rowvec<T>& a,const rowvec<T>& b){
		if(a.size() != b.size())
			throw "Illegal Inner Multiply!";
		T res = 0;
		int n = a.size();
		for(int i = 0; i < n; ++i) res += a[i] * b[i];
		return res;
	}
	friend T operator * <> (const rowvec<T>&, const colvec<T>&);
};

template <class T> class colvec : public matrix<T>{
public:
	inline colvec():matrix<T>(){}
	inline colvec(int _n):matrix<T>(_n,1){}
	inline T& operator [] (const int& i){
		return matrix <T> ::at(i);
	}
	inline T operator [] (const int& i)const{
		return matrix <T> ::at(i);
	}
	inline int size()const{
		return matrix <T> ::size().first;
	}
	friend T operator * (const colvec<T>& a, const colvec<T>& b){
		if(a.size() != b.size())
			throw "Illegal Inner Multiply!";
		T res = 0;
		int n = a.size();
		for(int i = 0; i < n; ++i) res += a[i] * b[i];
		return res;
	}
	friend T operator * <>(const rowvec<T>&, const colvec<T>&);
};

template <class T> inline T operator * (const rowvec<T>& a,const colvec<T>& b){
	if(a.size() != b.size())
		throw "Illegal Multiply!";
	T res=0;
	int n = a.size();
	for(int i = 0; i < n; ++i) res += a[i] * b[i];
	return res;
}

#endif
